#include <iostream>
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <string.h>
#include <stdio.h>
#include <fstream>

using namespace std;

unsigned int reducir(unsigned int n, double fred, unsigned int * img);
void ampliar(unsigned int n, unsigned int * img, unsigned char * imgCpy, unsigned int dim);
unsigned int indiceEnA(unsigned int fila, unsigned int columna, unsigned int* filaPtr, int n);
double sum(double* A, double* x, unsigned int desde, unsigned int hasta, unsigned int* colInd);

void SaltearComentarios(FILE *fp)
{
	int ch;
	char line[100];

	while ((ch = fgetc(fp)) != EOF && isspace(ch))
	{
	    ;
	}
	if (ch == '#')
	{
		fgets(line, sizeof(line), fp);
		SaltearComentarios(fp);
	}
	else
	{
        fseek(fp, -1, SEEK_CUR);
	}
}

int main(int argc, char *argv[])
{
    char* imgOrig = argv[1];
    unsigned int fred = 1;
    unsigned int reduccion = atoi(argv[2]);
    while(reduccion > 0)
    {
        fred = fred * 2;
        reduccion--;
    }
    double lambda = atof(argv[3]);

    char imgCopy[strlen(imgOrig) + 1];
    strcpy(imgCopy, imgOrig);


    // Lectura del PGM:
    FILE* fp = fopen(strcat(imgOrig,".pgm"),"rb");
    unsigned int height, width;
	int primerC, segC, max, c1, c2, c3, f;

    primerC = fgetc(fp);
	segC = fgetc(fp);
	SaltearComentarios(fp);
	c1 = fscanf(fp, "%d", &width);
	SaltearComentarios(fp);
	c2 = fscanf(fp, "%d", &height);
	SaltearComentarios(fp);
	c3 = fscanf(fp, "%d", &max);
	cout << "Altura:" << height << endl << "Ancho:" << width << endl;
    unsigned int* image = new unsigned int[height*width];

	if (primerC != 80 || segC != 53 || c1 != 1 || c2 != 1 || c3 != 1 || max > 255)
    {
        cout << c3 << " " << max << " " << "La imagen seleccionada no tiene formato PGM.\n" << endl;
        return 1;
    }

	fgetc(fp);//Paso de largo el caracter siguiente ya que no me interesa (segun especificacion, generalmente es un newline)
    unsigned char temp;
	for (int i = 0; i < height*width; i++)
	{
            temp = fgetc(fp);
            image[i] =  temp;
	}

	fclose(fp);


	//Reduccion:

    unsigned int dim = reducir(height, fred, image);


    //Escritura PGM:
	FILE* fpOut = fopen(strcat(imgCopy,"OutR.pgm"),"wb");

    fprintf(fpOut, "P5 ");
    fprintf(fpOut, "%d %d ",dim,dim);
    fprintf(fpOut,"%d ",max);

    unsigned char imageCpy[dim*dim]; // aca guardaremos la imagen filtrada (sin ampliar)

    for (f = 0; f < dim*dim; ++f)
    {
        fputc(image[f], fpOut);
    }

    fclose(fpOut);


    //Filtrado:

    //matriz del sistema K * x = b
    int n = dim; // dim = cant filas/columnas de la imagen 'image' (reducida)

    //vector columna independiente (b)
    //copiamos la imagen en b (b contiene las condiciones iniciales)
    int dimb = n*n;
    double* b = new double[dimb];
    for(int i=0; i<dimb; i++){
        b[i] = image[i];
    }

    //Creamos la matriz A (A contiene la banda de la matriz, o sea todo sin los triangulos sup e inf de ceros, junto con el vector b al costado derecho)
    int triangCeros = 0;
    for(int i=1; i<=n*n-n-1; i++){
        triangCeros += i;
    }
    int dimA = (n*n)*(n*n) - 2*triangCeros + dimb;
    double* A = new double[dimA];
    for(int i=0; i<dimA; i++) {A[i] = 0;}

    //pongo el vector b
    int ib = 0; //indice para recorrer b
    int j = 1;
    int k = n+3;
    for(int count=1; count<=n+1; count++){
        A[n+j] = b[ib];
        ib++;
        j+=k;
        k++;
    }
    int alturaBloquesTotal = (n-2)*(n-2) + 2*(n-3); //altura de toda la banda excepto los 2 triangulos de unos
    int index = n+j-1; // el valor de b siguiente se encuentra justo despues de la primera fila de los 5NoNulos
    for(int count=1; count<=alturaBloquesTotal; count++){
        A[index] = b[ib];
        index += (2*n)+2;
        ib++;
    }
    k = (2*n)+1;
    for(int count=1; count<=n+1; count++){
        A[index] = b[ib];
        ib++;
        index += k;
        k--;
    }

    //pongo el triangulo superior izquierdo de 1's
    int i0 = 0;
    A[i0] = 1;
    for(int i=(n+2)+1; i<=(n+2)+n; i++){
        A[i0+i] = 1;
        i0+=i;
    }

    //pongo el triangulo inferior derecho de 1's
    i0 = dimA-2;
    A[i0] = 1;
    for(int i=(n+2)+1; i<=(n+2)+n; i++){
        A[i0-i] = 1;
        i0-=i;
    }

    //primero algunos offsets:
    unsigned int inicio5NoNulos = 0; //donde empieza la 1era fila de -1...-1 5 -1...-1
    unsigned int offset5NoNulos; //lo que ocupa un bloque de renglones 5NoNulos
    unsigned int inicioBloque2Unos; //donde empieza la 1era fila que contiene el bloque de dos 1's
    unsigned int offsetBloque2Unos; //lo que ocupa un bloque de renglones 2Unos
    //ahora los inicializamos bien:
    for(int i=1; i<=n+1; i++) {inicio5NoNulos+=(n+i)+1;} //el +1 es por el valor del b_i
    offset5NoNulos = (n-2)*(2*n+2);
    inicioBloque2Unos = inicio5NoNulos + offset5NoNulos;
    offsetBloque2Unos = 2*(2*n+2);

    //pongo los Bloques2Unos
    index = inicioBloque2Unos;
    for(int count=1; count<=n-3; count++){
        A[index + n] = 1;
        A[index + (2*n+2) + n] = 1;
        index += offsetBloque2Unos + offset5NoNulos; // llego al prox bloque2Unos
    }

    //pongo los Bloques5NoNulos
    index = inicio5NoNulos;
    for(int count=1; count<=n-2; count++){
        for(int count2 =1; count2<=n-2; count2++){
            A[index] = -1/lambda; //arriba
            A[index + (n-1)] = -1/lambda; //izquierda
            A[index + (n-1) + 1] = (lambda+4)/lambda; // centro
            A[index + (n-1) + 1 + 1] = -1/lambda; // derecha
            A[index + (n-1) + 1 + 1 + (n-1)] = -1/lambda; // abajo
            index += 2*n+2;
        }
        index += offsetBloque2Unos; // llego al prox bloque5NoNulos
    }


    //Triangulamos la matriz A

    /* Primero creamos dos vectores: col_ind y fila_ptr.
     * Sea K la matriz banda del sistema completa (de dimension dim^4).
     * col_ind[i] es la columna en K del elemento A[i].
     * fila_ptr[j] me dice que en A[col_ind[j]] empieza una nueva fila de K.
     * col_ind tiene la misma dimension que A.
     * fila_ptr tiene dimension igual a dim^2+1.
     * Por convencion, la ultima posicion de fila_ptr contiene (dimA+1).
     * (Guardamos el tamanio de la banda de K en la ultima posicion de fila_ptr).
     */

     int dimCol_Ind = dimA;
     unsigned int* col_ind = new unsigned int[dimCol_Ind];
     int dimFila_Ptr = (n*n)+1;
     unsigned int* fila_ptr = new unsigned int[dimFila_Ptr];

     //relleno col_ind
     //pongo las columnas correspondientes al triangulo superior izquierdo de 1's
     int c = 0; //indice para recorrer col_ind
     j = 0; //a cada elemento de cada fila del triangulo superior izquierdo de 1's le corresponde la secuencia de columnas en K de 0 a n+j
     for(int count=1; count<=n+1; count++){
        for(int i=0; i<=n+j; i++){
            col_ind[c] = i;
            c++;
        }
        j++;
        col_ind[c] = n*n; //todos los b_i estan en la columna n*n de la matriz K
        c++;
     }
     //pongo las columnas correspondientes al grupo de bloques 5NoNulos y 2Unos
     int kl = 1; //a cada elemento de cada fila del grupo de bloques 5NoNulos y 2Unos le corresponde la secuencia de columnas en K de kl a 2n+kr
     int kr = 1;
     for(int count=1; count<=alturaBloquesTotal; count++){
        for(int i=kl; i<=(2*n)+kr; i++){
            col_ind[c] = i;
            c++;
        }
        kl++;
        kr++;
        col_ind[c] = n*n;
        c++;
    }
    //pongo las columnas correspondientes al triangulo inferior derecho de 1's
     j = 1; //a cada elemento de cada fila del triangulo inferior derecho de 1's le corresponde la secuencia de columnas en K de alturaBloquesTotal+j a 2n+alturaBloquesTotal+1
     for(int count=1; count<=n+1; count++){
        for(int i=alturaBloquesTotal+j; i<=(2*n)+alturaBloquesTotal+1; i++){
            col_ind[c] = i;
            c++;
        }
        j++;
        col_ind[c] = n*n;
        c++;
     }

    //relleno fila_ptr
    //pongo las filas respecto del triangulo superior izquierdo de 1's
    f = 0; //indice para recorrer fila_ptr
    j = 0;
    for(int count=1; count<=n+1; count++){
        fila_ptr[f] = j;
        f++;
        j+=n+count+1; //el +1 es por los b_i
    }
    //pongo las filas respecto del grupo de bloques 5NoNulos y 2Unos
    for(int count=1; count<=alturaBloquesTotal; count++){
        fila_ptr[f] = j;
        f++;
        j+=(2*n)+2; //2n+2 es la cantidad de elementos en la fila que tiene 5 no nulos o la del 1 del medio (+ el b_i)
    }
    //pongo las filas respecto del triangulo inferior derecho de 1's
    for(int count=n+1; count>=1; count--){
        fila_ptr[f] = j;
        f++;
        j+=n+count+1; //el +1 es por los b_i
    }
    //guardo por defecto al final el tamanio de la banda de K (la dimension de A)
    fila_ptr[dimFila_Ptr-1] = dimA;

    //Por comodidad para cuando triangulemos, nos guardamos los indices en A de la diagonal principal
    int dimDiagonal = n*n;
    double* diagonal = new double[dimDiagonal];
    int d = 0; //para recorrer diagonal
    int iA = 0; //indice en A
    for(int count=1; count<=n; count++){
        diagonal[d] = iA;
        d++;
        iA += (n+2)+count;
    }
    for(int count=1; count<=alturaBloquesTotal+1; count++){
        diagonal[d] = iA;
        d++;
        iA += (2*n)+2;
    }
    for(int count=n; count>=0; count--){
        diagonal[d] = iA;
        d++;
        iA += (n+2)+count;
    }

    //Tambien nos guardamos los indices en A de los valores b_i en orden
    int dimBPos = n*n;
    unsigned int* bPos = new unsigned int[dimBPos];
    for(int i=1; i<dimFila_Ptr; i++){
        bPos[i-1] = fila_ptr[i]-1;
    }

    //muestra A antes de triangular:
    /*cout << endl;
    cout << "A antes de triangular: ";
    for(int j=0; j<dimA; j++){
        cout << A[j] << " ";
    }*/

    //muestra valores de b_i antes de la triangulacion
    /*cout << endl;
    cout << "b_i antes de triangular: ";
    for(int j=1; j<dimFila_Ptr; j++){
        cout << A[ fila_ptr[j]-1 ] << " ";
    }*/


    //Triangulamos A:
    /* Por como es la matriz A, sabemos que todos los elementos de su diagonal son no nulos,
     * luego no hay necesidad de intercambio de filas. */
    double ai; //aca guardamos los nuevos valores de los A_i
    int p = 0;
    unsigned int pNumM; //POSICION en A del numerador del multiplicador (el primer elem de la fila q estoy viendo)
    double numM; //numerador del multiplicador
    unsigned int indexPivote; //POSICION en A del denominador del multiplicador (el pivote actual)
    d = 0; //para recorrer el vector diagonal
    int fdesde = 1; //fila desde (en cada iteracion de gauss)
    int fhasta = n+1; //fila hasta (en cada iteracion de gauss)
    int cdesde; //col desde (en cada iteracion de gauss)
    int chasta = n; //col hasta (en cada iteracion de gauss)

    for(cdesde=0; cdesde<n*n; cdesde++){ //pasos de la eliminacion gaussiana (la columna n*n corresponde a los b_i)
        indexPivote = diagonal[d]; //indice en A del elemento de la diagonal del paso 'col' de gauss
        int indexbPos = 1+p; //indice en A del b_i inmediatamente debajo del b_i de la fila pivote
        for(int fila=fdesde; fila<fhasta; fila++){ //trabajo con cada fila
            int indexFilaPasoCol = indexPivote; //elementos de la fila de paso 'col' de gauss
            pNumM = indiceEnA(fila,cdesde,fila_ptr,n); //me guardo el numerador del multiplicador para esta fila "fila"
            numM = A[pNumM]; //porque despues piso A[pNum]
            for(int col = cdesde; col<=chasta; col++){
                int i = indiceEnA(fila,col,fila_ptr,n);//posicion del elemento a modificar
                ai = A[i] - ((numM / A[indexPivote]) * A[indexFilaPasoCol]); //comb lineal de renglones
                A[i] = ai;
                indexFilaPasoCol++;
            }
            //resto por separado los valores de b_i
            unsigned int bpi = bPos[indexbPos];
            ai = A[bpi] - ((numM / A[indexPivote]) * A[indexFilaPasoCol]);
            A[bpi] = ai;
            //incrementos
            indexbPos++;
        }

        fdesde++; //voy al siguiente paso de gauss
        if (cdesde < (n*n)-n-1){ //cuando llego al triangulo inferior derecho de 1's dejo de hacer fhasta++ y chasta++
            fhasta++;
            chasta++;
        }
        d++; //avanzo al prox elemento de la diagonal
        p++; //avanzo al prox elemento del vector b
        //k++; //paso a la posicion del elemento inmediatamente debajo del pivote
    }


    //Encontramos la solucion (sustitucion hacia atras)
    int dimX = n*n;
    double* x = new double[dimX]; //vector solucion (pixeles filtrados de arriba abajo, izquierda derecha)

    for(int i=0; i<dimX; i++) {x[i] = 0;}

    for(int i=dimX-1; i>=0; i--){
        double bi = A[ bPos[i] ];
        double aii =  A[ indiceEnA(i,i,fila_ptr,n) ];
        unsigned int desde = indiceEnA(i,i+1,fila_ptr,n); // posicion en A del termino a_{i,i+1} . x_{i+1}
        unsigned int hasta;
        if(i==0){
            hasta = 0; //para evitar que 'hasta' me de negativo
        }else{
            hasta = fila_ptr[i]-2; // posicion en A del termino a_{i,dimX-1} . x_{dimX-1}
        }
        x[i] = (bi - sum(A,x,desde,hasta,col_ind)) / aii;
    }

    //muestra A luego de triangular:
   /* cout << endl;
    cout << "A luego de triangular: ";
    for(int j=0; j<dimA; j++){
        cout << A[j] << " ";
    }*/

    //muestra x (la solucion obtenida):
    /*cout << endl;
    cout << "vector solucion (x): ";
    int fueraRango = 0;
    for(int j=0; j<dimX; j++){
        cout << x[j] << " ";
        if (x[j]<0 || x[j]>255) {fueraRango++;}
    }*/
   // cout << "cant de pixeles fuera de rango: " << fueraRango << endl; //siempre da 0! :) O sea q esta OK!

    //muestra bPos
    /*cout << endl;
    cout << "bPos: ";
    for(int j=0; j<dimBPos; j++){
        cout << bPos[j] << " ";
    }*/

    //muestra col_ind
   /* cout << endl;
    cout << "col_ind: ";
    for(int j=0; j<dimCol_Ind; j++){
        cout << col_ind[j] << " ";
    }*/

    //muestra fila_ptr
   /* cout << endl;
    cout << "fila_ptr: ";
    for(int j=0; j<dimFila_Ptr; j++){
        cout << fila_ptr[j] << " ";
    }*/

    //muestra posiciones de los elementos en la diagonal en A:
    /*cout << endl;
    cout << "diagonal: ";
    for(int j=0; j<dimDiagonal; j++){
        cout << diagonal[j] << " ";
    }
    */
    //cout << "  porcentaje elementos en banda: " <<  ((double) dimA / (dim^4) * 100) << "   ";

    //muestra valores de b_i luego de la triangulacion
   /* cout << endl;
    cout << "b_i luego de triangular: ";
    for(int j=1; j<dimFila_Ptr; j++){
        cout << A[ fila_ptr[j]-1 ] << " ";
    }*/
//---------------------------ACA TERMINA EL FILTRADO--------------------------------------------


    //Pasamos el vector solucion (x) (todavia sin ampliar) a 'image'
    for(int i=0; i<dimX; i++){
        image[i] = x[i];
    }


    //Ampliado:

    for(int i = 0; i < dim*dim; i++)
    {
        unsigned char temp = (unsigned char) image[i];
        imageCpy[i] = temp;
    }

    ampliar(height,image,imageCpy,dim);


    //Escritura PGM:
	FILE* fpOut2 = fopen(strcat(imgCopy,"Out.pgm"),"wb");

    fprintf(fpOut2, "P5 ");
    fprintf(fpOut2, "%d %d ",height,height);
    fprintf(fpOut2,"%d ",max);

    for (f = 0; f < height*height; ++f)
    {
        fputc(image[f], fpOut2);
    }

    fclose(fpOut2);


    //Libero memoria dinamica:
    delete[] x;
    delete[] bPos;
    delete[] diagonal;
    delete[] fila_ptr;
    delete[] col_ind;
    delete[] A;
    delete[] b;
    delete[] image;

    return 0;
}

unsigned int reducir(unsigned int n, double fred, unsigned int * img)
{

    unsigned int dim = n/fred;//dimension de la imagen reducida

    unsigned int k = floor(n/dim);//filas o columnas a borrar

    unsigned int c=0,f=0;

    for(int ft=0; ft < dim; ft++)//recorre filas de T
    {
        for(int ct=0; ct < dim; ct++)//recorre columnas de T
        {
            img[ft*dim+ct] = img[f*n+c];//llena T con valores de imagen
            if((c+k) < n)
            {
                c = c + k;
            }
        }
        if((f+k)<n)
        {
            f = f + k;
        }
        c = 0;
    }

    return dim;
}

void ampliar(unsigned int n, unsigned int * img, unsigned char * imgCpy, unsigned int dim)
{

    unsigned int k = floor(n/dim);//filas o columnas a borrar

    unsigned int c=0,f=0;

    for(int ft=0; ft < n; ++ft)//recorre filas de T
    {
        for(int ct=0; ct < n; ++ct)//recorre columnas de T
        {
            unsigned char temp = imgCpy[f*dim+c];
            img[ft*n+ct] = temp;//llena T con valores de imagen
            for(int r = 0;r<k-1;r++)
            {
                ct++;
                temp = imgCpy[f*dim+c];
                img[ft*n+ct] = temp;
            }

            c++;

        }
        //COPIO LA FILA ANTERIOR
        unsigned int ftAnt = ft;

        for(int r = 0;r<k-1;r++)
        {
            ft++;
            for(int cj=0;cj < n; cj++)
            {
                    unsigned char temp = img[ftAnt*n+cj];
                    img[ft*n+cj] = temp;
            }
        }


        f++;

        c = 0;
    }

}

unsigned int indiceEnA(unsigned int fila, unsigned int columna, unsigned int* filaPtr, int n){
//dada la fila "fila" y la columna "columna", me devuelve el indice del elemento en A que
//se encuentra en esas coordenadas
    int offset;
    if(fila<=n){
        offset = 0;
    }else{
        offset = fila-n;
    }
    return (filaPtr[fila] + columna - offset);
}

double sum(double* A, double* x, unsigned int desde, unsigned int hasta, unsigned int* colInd){
//calcula la suma de los terminos despues del pivote
    double sumatoria = 0;
    while(desde<=hasta){
        double termino = A[desde] * x[ colInd[desde] ];
        sumatoria += termino;
        desde++;
    }
    return sumatoria;
}










